rm(list=ls()) 
library(data.table)
library(stargazer)
library(lfe)
library(dplyr)
library(tidyr)
library(haven)
# Make all the tables 
# Read in fles
setwd("C:/Users/mg5797/Documents/Pollution/ReplicationPackage")

# import temperature: 
tempdf=read.csv("Derived Data/weatherdf.csv",header=TRUE,stringsAsFactors = FALSE)

tempvarsavg="ay_mean_snow  +aymean_high_num_snow+aymean_severe_num_snow+ay_avg_mean_max + ay_avg_mean_min  + ay_avg_bint80  + ay_avg_bint90  + ay_avg_bint100  + ay_avg_bint20 + ay_avg_bint10+ ay_avg_bint0+ ay_mean_prcp +aymean_high_num_prcp+aymean_severe_num_prcp"

districtcovariates="pct_hisp_2001:factor(year) + pct_black_2001:factor(year) + pct_male_2001:factor(year) + pct_white_2001:factor(year) + enroll_2001:factor(year) + pct_ell_2001:factor(year) + pct_speced_2001:factor(year) + pct_FRPM_2001:factor(year) + missing_pct_ell_2001 + missing_pct_speced_2001 + missing_pct_FRPM_2001 "

acs="pctemploy2001:factor(year)  + pctsinglemother2001:factor(year) + pctmanufac2001:factor(year) + pctutility2001:factor(year) + pctbachelorhigher2001:factor(year)"




# loop through all the distances 

output_tables <- function(dfbartik, distance,monthfilter,subjectfilter,het){
  
  filetoimport=paste("Derived Data/dfbartik2001_",distance,"kmanalysis.csv",sep="")
  
  dfbartik=fread(filetoimport)
  
  dfbartik=merge(dfbartik,tempdf,by=c("leaid","year"),all.x=TRUE)
  
  dfbartik$pctemploy2001=as.numeric(dfbartik$pctemploy2001); dfbartik$pctsinglemother2001=as.numeric(dfbartik$pctsinglemother2001); dfbartik$pctmanufac2001=as.numeric(dfbartik$pctmanufac2001); dfbartik$pctutility2001=as.numeric(dfbartik$pctutility2001)

  dfbartik$pctbachelorhigher2001=as.numeric(dfbartik$pctbachelorhigher2001)
    ############### Bartik 2001 IV ################# 
  if(monthfilter=="May"){
    dfbartik_may=dfbartik[dfbartik$testmonth=="May",]
    
  }
  else{
    dfbartik_may=dfbartik
    
    
  }
  
  
  if(subjectfilter=="rla"){
    dfbartik_may=dfbartik_may[dfbartik_may$subject=="rla",]
    
    
  }
  else if(subjectfilter=="mth"){
    dfbartik_may=dfbartik_may[dfbartik_may$subject=="mth",]
    
  }
  else{
    dfbartik_may=dfbartik_may
  }
  
  dfbartik_may=dfbartik_may %>% drop_na(cs_mn_all,avgpm25_9, cs_mn_all_lag_cohort) # same drop conditons as last sample
  nrow(dfbartik_may)
  

  # drop missing
  dfbartik_may=dfbartik_may %>% drop_na(cs_mn_all,avgpm25_9, cs_mn_all_lag_cohort, pct_hisp_2001,pct_male_2001,enroll_2001,pct_ell_2001,pct_speced_2001,pct_FRPM_2001)
  
  nrow(dfbartik_may)
  
  # instrument: 
  instname=paste("bartikleave_",distance,sep="")
  # weights:
  weightval=paste("totalprod_fuel",distance,sep="")
  
  if(het==0){
    
    
    bartik_none=as.formula(paste("cs_mn_all~ cs_mn_all_lag_cohort + ", districtcovariates ,"  |factor(year) +factor(subject) + factor(grade)+factor(leaid) |(avgpm25_9 ~ ",instname,")|leaid",sep="" ))
    
    resultbartik_none=felm(bartik_none,dfbartik_may,weights=dfbartik_may[[weightval]],na.action=na.omit)
    
  #  district covariates
  bartik_dist=as.formula(paste("cs_mn_all~ cs_mn_all_lag_cohort + ", districtcovariates ,"  |factor(year) +factor(subject) + factor(grade)+factor(leaid) |(avgpm25_9 ~ ",instname,")|leaid",sep="" ))
  
  resultbartik_dist=felm(bartik_dist,dfbartik_may,weights=dfbartik_may[[weightval]],na.action=na.omit)

  
  bartik_sort=as.formula(paste("cs_mn_all~ cs_mn_all_lag_cohort + ",districtcovariates, "+ ", acs ," |factor(year) +factor(subject) + factor(grade)+factor(leaid) |(avgpm25_9 ~",instname,")|leaid",sep="" ))
  resultbartik_sort=felm(bartik_sort,dfbartik_may,weights=dfbartik_may[[weightval]],na.action=na.omit)

  
  # district + acs +weather long 
  bartik_weathavg=as.formula(paste("cs_mn_all~ cs_mn_all_lag_cohort + ", districtcovariates, " + ", acs ," + ", tempvarsavg ,"  |factor(year) +factor(subject) + factor(grade)+factor(leaid) |(avgpm25_9 ~ ",instname,")|leaid",sep="" ))
  resultbartik_weathavg=felm(bartik_weathavg,dfbartik_may,weights=dfbartik_may[[weightval]],na.action=na.omit)
 
  # export 
  fileout= paste("output/Table_bartik2001leave_",distance,"km_month",monthfilter,".tex",sep="")
  
  stargazer(resultbartik_dist,resultbartik_sort,resultbartik_weathavg, keep="avgpm25_9",digits=4,
               out=fileout )

}
  
  
  
}




output_tables(dfbartik,40,"May","all",0)

